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We propose the consistent statistical approach to consider a wide class of classical open systems 
whose states are specified by a set of positive integers(occupation numbers). Such systems are often 
encountered in physics, chemistry, ecology, economics and other sciences. Our statistical method 
based on ideas of quantum theory of open systems takes into account both discreteness of the 
system variables and their time fluctuations - two effects which are ignored in usual mean field 
dynamical approach. The method let one to calculate the distribution function and (or) all moments 
of the system of interest at any instant. As descriptive examples illustrating the effectiveness of the 
method we consider some simple models:one relating to nonlinear mechanics, and others taken from 
population biology .In all this examples the results obtained by the method for large occupation 
numbers coincide with results of purely dynamical approach but for small numbers interesting 
differences and new effects arise. The possible observable effects connected with discreteness and 
fluctuations in such systems are discussed. 
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I. INTRODUCTION 

Among the vast number of classical open dynamical 
systems under consideration in physics, chemistry, biol- 
ogy, economics and other sciences there are many such 
whose states in accordance with the sense of the problem 
are specified by a set of integer variables {{rii}) , where i 
= 0,1,2,..N (N-number of degrees of freedom). For exam- 
ple in physics n^- are occupation numbers of cell states 
in phase space, in chemistry - numbers of molecules of 
reactive elements, in ecology -numbers of individuals in 
populations which live in the area and interact with other 
populations, in economics the number of companies oper- 
ating on the market. Since all these systems are classical 
their dynamics as a rule is described by a system of dif- 
ferential equations of the form 






F^ ({n„}) 



(1) 



where Fi{{no}) -some nonlinear functions, depending on 
concrete problem. Obviously notation Eq. ([ij implies 
that variables n^ in system Eq. ^ -are considered as 
continuous. In the case when all n^ ^ 1 such approach 
can be easy justified. From the physical point of view 
the system of equations ([T]) corresponds to mean field 
approximation and ni are occupation numbers averaged 
over some appropriate statistical ensemble. In the case 
when all ?ii ^ 1, dynamical description becomes inade- 
quate and the question naturally arises: is there consis- 
tent statistical approach which takes into account both 
discreteness of variables Ui and their time fluctuations 
that may be not small. In addition it is naturally to de- 
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mand that such approach gave the same results as dy- 
namical description in the large rii limit. In this paper 
we propose such approach based on the ideas of quan- 
tum theory of open systems (QTOS) and consider some 
examples that demonstrate its effectiveness. The rest of 
the article organized as follows. In the Sect. 2 we briefly 
describe minimal information from (QTOS) which is nec- 
essary for understanding of the method used and present 
main steps of our method. In Sect. 3 we consider simple 
model of nonlinear autonomous oscillator with soft excit- 
ing mode and give its statistical description on the basis 
of the method proposed. All the main features of the 
method clearly come to light already in this representa- 
tive example. In Sect. 4 with the help of our approach 
we consider some problems from population dynamics 
relating to evolution of two interaction populations liv- 
ing in a certain area.. We show that in the case of small 
populations statistical description leads to a number of 
differences from the ordinary dynamical picture. On the 
other hand in the case of large occupation numbers both 
descriptions are virtually identical. In conclusion we dis- 
cuss some generalizations of the method and its possible 
experimental veriflcation. 



II. DESCRIPTION OF THE METHOD 

In this section we briefly remind the main points of the 
method proposed by author earlier [l| which allows one 
to make the transition from known dynamical equations 
of classical open system to the master equation for its 
quantum analogue. The method based on the correspon- 
dence that can be set between quantum master equation 
in the Lindblad form and the Liouville equation for distri- 
bution function in phase space of the classical system of 
interest. This correspondence allows one using classical 



equations of motion to restore the form of all operators 
involved in the Lindblad equation. Thereby we can ap- 
ply the procedure of quantization at least in semiclassical 
approximation in the case of a large class of nonhamil- 
tonian dynamical systems. In short (all details see in 
[i| ) the recipe of quantization proposed consists of three 
consecutive steps. 

Stepl: The input classical dynamical equations 
should be presented in the form allowed the quantiza- 
tion (FAQ). For the purposes of present paper the most 
convenient form is complex representation of equations 
of motion: 



dzj 
Itt 



dH 



EZ-^dRa dR, 

I -Kq, —r-:_ R 



dz* 



dz* 



(2) 



where 



V2 



V2 



are complex dynami- 



cal coordinates of the system of interest, H, Ra, Ra are 
functions of Zi, z* ( H is real function, and Ra, Ra are 
complex, Ra means function which conjugate to Ra)- It 
is necessary to emphasize that it is the most delicate step 
of the method because it is difficult exactly to formalize 
this point. 

Step2. Having in hands representation Eq. ([2]) we 
can use classical function H, Ra, Ra and with their 
help determine their quantum analogues - operators H, 
Ra, i?+. For this purpose the variables Zi, z* must 
be replaced by the correspondence Bose operators a^ 



and af with usual commutation rules: [ai,a^] = 6 



ij, 



[ai,a.j 



= 0. 



Step3. The operators H, Ra, Ra found in this man- 
ner should be substituted into the quantum Lindblad 
equation for the evolution of density matrix of the sys- 
tem: 



dp 
'dt 



i?,p] +J2[Ra-?,R^ 



Ra , p ■ Ro 



The correspondence principle guarantees us that such ap- 
proach will give correct description of evolution of quan- 
tum open system at least with the accuracy up to the 
first order in h. 



III. REPRESENTATIVE MODEL. 

AUTONOMOUS NONLINEAR OSCILLATOR 

WITH SELF EXCITING MODE 

It is convenient to demonstrate on concrete example 
all the features of the approach proposed. We consider 
an oscillator with nonlinear damping in situation when 
its equilibrium point loses stability and small fiuctua- 
tions switch the system to the new stationary state that 
corresponds to the closed trajectory (limit cycle). The 
following system of equations gives the correct mathe- 
matical description of the behavior of the oscillator in 
the vicinity of bifurcation point (see [J]): 



dx 

dy 
dt 



ujy + px — X (x^ -I- y^) , 
—ojx + py — y (x + y ) 



(4) 



where x, y are coordinates of oscillator in phase space 
and w -its frequency. The system of the equations of 
motion Eq. Q can be written in the complex form as 
one equation 



dz 



-iujz + pz — 2z \z\ 



(5) 



where z = ^^p . It is easy to show that equation Eq. 

([S]) can be represented in the FAC. For this purpose we 
introduce the functions H = ujz* z, Ri = ^/pz* and R2 = 
z^ . One can verify by direct checking that r.h.s. of Eq. 
([5]) can be rewritten in the form 



dz 
Iti 



dH 

oz* 



Ri 



dRi 

dz* 



Ri 



dRi 
dz* 



R: 



■dR2 
dz* 



R 



dR2 
dz* 

(6) 
According to the recipe of quantization the Lindblad 
equation for quantum analog of the system Eq. Q can 
be written in the form 



dp 
'di 



= -i [H, p] + [Rip, R+] + [R2P, R+] + h. 



c, 



(7) 



where H = Lod^a, Ri = ^/pd^ and R2 —d^. 

We are interesting in the stationary solution of Eq. [7] 
and from physical considerations we will seek it in the 
form: p — 'p{N) — '^ \n) p^^ {n\ , where \n) are the eigen- 
functions of the operator N = d^d. It is convenient to in- 
troduce the generating function G{u) for the coefficients 



(3) p„. By definition G{u) 



Y^ p„ • u".One can obtain 

n=0 

from Eq. [7] in stationary case the following equation for 
the function G{u): 



(l + u) 



fG 

dv? 



dG 
pu—- - pG{u) = 0. (8) 

du 



The solution of Eq. |S] that satisfies all conditions of the 
problem may be represented as 



G{u) 



<^{l,p,p{l + u)) 
Hl,p,2p) - 



(9) 



where ^{a,c,x) is well known confiuent hypergeometric 
function (see [|) - ^a,c,x) = 1 + fi + ^^x^ -h .... 

Note that condition G{u) = 1 corresponds to normaliza- 
tion of p^^{n)— namely ^p„ = 1. Having in hands the 

n 

expression Eq. [9] for the generation function one can find 
the average values for any physical quantity of interest 
in the stationary state that is all moments of the distri- 
bution p(n). For example moments of first and second 



order are determined by relations: 



E " ■ Pn 



dG 



d^G 



n^ — n = ^^- (all derivatives are taken at point u = 1). 
Let us consider now the behavior of our system in two 
limiting cases: /i ^ 1 and /i <C 1. When u ^ 1 using 
the asymptotic formula for $(a, c, a:) (see 0)we obtain 
for G (u) the next expression 



Gin) 



M^-'^) 



1 



(10) 



The average number of quanta in stationary state n and 
their dispersion a = n'^ — n^ in this case are: n — ^^-^ 
and (T = ^'^^^ ■ The relative fluctuation of quanta gen- 
erated in stationary state is -^ = ^^^^ tends to zero 
when /i ^ 1. The obtained result testifies validity of 
deterministic approach in this case since in classical case 
Uci — \z\ = ^ when the system is moving along the limit 
cycle. Now let us consider the opposite case when ^ ^ 1. 
Using the series expansion for $(a, c, x) 0] we obtained 
for G{u) the next expression (up to the first order in fi) 



Giu) 



(2 + u){l + fi) + ii{l + uf 
3 + 7Ai 



(11) 



The relation Eq. [TT] implies that G{u) tends to Go{u) = 
^^ when fi tends to zero. The first two moments of 
the distribution p(n) are: n = ^ and ct = | . The relative 

fluctuation of n = -^ is equal to \/2~ in this case. We see 
that in contrast to classical situation n ^ when /i tends 
to zero moreover fluctuations of occupation number turn 
out to be large and very essential. 

It is necessary to note that in classical system Eq. |3] 
in addition to variable | z \ we have phase variable ip that 



satisfies to equation 



dip 
It 



However phase dynam- 



ics determines only velocity along limit cycle but does 
not influence on stationary state itself. In principle we 
can turn w to zero (which implies that H = 0) and all 
foregoing results do not change. It is the point which 
explains why the Lindblad equation can be successfully 
applied for statistical description of classical open sys- 
tems with integer variables. The reason is that in the 
case when H — in Eq. SI h may be eliminated from the 
Lindblad equation and wave properties of quantum sys- 
tem turn out to be irrelevant. After this crucial remark 
we can apply our approach to different classical systems 
with integer variables in particular population dynamics 
models. 



IV. STATISTICAL DESCRIPTION OF 
POPULATION DYNAMICS MODELS. 

We begin our consideration with well known Lotka- 
Volterra model (LVM) (see [4i]) describing "the interac- 
tion" between two populations: prays and predators (for 
example hares and lynxes) living in the same territory. 



Let ni(i) and n2(i) the current numbers of preys and 
predators correspondently. Then input dynamical equa- 
tions of the LVM can be written as 



— — = ni{a~ 2^2) , 
at 

— = -nAb-2n,)., 



(12) 



where coefficients a, b are positive and have clear ecolog- 
ical meaning. Note that for convenience we have chosen 
time unit thus to put coefficient of the term nin2 in Eq. 
[T^ equal to 2. Now let us demonstrate that system Eq. 
[12] can be represented in FAQ. For this purpose we intro- 
duce two auxiliary complex variables zi and Z2 so that 
ni = \zi\ and 712 = \z2\ ■ Let us assume that evolu- 
tion zi , Z2 in time is governed by the following system of 
equations: 



— = \Zl - Zl \Z2\ ; 



(13) 



dZ2 

~dt 



-\lz2 + Z2\zi\' 



It is easy to see that system Eq. [T3] implies the following 
equations for ni, n2 



drii 2 

—— = 2X^111 - 2^1712, 

dn2 2 

—— = -2A2n2 -I- 2nin2, 



(14) 



which is exactly coincides with Eq. [T2]if one lets a = 2A]^, 
b = 2\l. 

Now if we introduce three functions i?i — Xiz^^ R2 — 
\2Z2^ and i?3 = z\Z2 it is easy to verify that Eq. 1131 can 
be represented in the form 



dz\ 
~di 

dZ2 

dt 






i=l 
3 



(15) 



V- ijr dR^ p dR, 



dz* 



dz* 



Using the foregoing recipe of quantization we maintain 
that description of LVM which takes into account the 
discreteness of variables rii and n2 and their time fluctu- 
ations is given by the next master equation 



dp 



3 

E 



[R^P.Rt 



(16) 



where i?i = AiaJ", R2 = A2a2 and R3 ~ aiOj . Again 
we will interesting in only the solutions of Eq. [!§] which 
have the form 'p — Yl l'^i'^2) Pn^ n^ ('^i''^2|- In this case 

for the coefficients of the expansion p„j„2 ^^ obtain the 
next general equation 



dPn 



dt 



= 2Ai [nip„^_i_„^ - (m + 1) p„^„J +2A2 (n2 + 1) p„^ 



,n2 + l "-aPnjnj 



r 



+2 [(m + 1) n2/?„^+i^„^_i - (^2 + 1) nip„i, 

OZ) 



It is convenient to introduce the generating function 



G{u,V,t) = J2 Pmn2 



.Then after the simple 



algebra we find that Eq. [T7] implies the next equation for 
G{u,v,t): 



-m = 2A?(z.-l)-(^G)+2A^(l-t;) — +2(x,-«)- 



92 



{vG). 



dudv 

(18) 

The Eq. [17] and Eq. [18] in principle give us all the nec- 
essary information about statistical behaviour of LVM. 
Now let us consider concrete results following from these 
equations. Note that even in stationary case it is difficult 
to find exact analytical solutions of Eq. [18] for arbitrary 
Aland A2. But in the special case when A2 = 1 + A^ such 
solution easily can be found and have the form G'(w, v) — 

where k = ttt2- Since G{u, v) — g{u)-g{v) 



(.i-xT 



it is clear that mand 712 are independent variables in 
this case and ni = ri2 — j^ — Aj^.The dispersion erf in 

this case is equal to A^ + Aj^ and if we calculate relative 
fluctuation of ni in stationary state we get the result: 

/5T 



Sni 



1 + 



Thereby we see this quantity 



is not small and can be easily measured. More detailed 
analysis of statistical behaviour of LVM following from 
Eq. [T7] for arbitrary Ai and A2 will be carried out else- 
where. Here we show only that using Eq. [18] one can 
easy obtain a collection of explicit relations connecting 
different moments of distribution p(ni, 712). For example 
if we differentiate stationary Eq. [T8]with respect to u and 
after that put u — v — 1 we obtain the simple relation 
between moments of first and second order which reads 



as: 



Xi (1 + ni) — ni — nin2 — 0. 



(19) 



In a similar manner by differentiating of Eq. [18] with 
respect to v we obtain the second relation: 



\2 — 

A2n2 



+ ni + nin-i — 



(20) 



Relations Eq. [TOland Eq. 
connecting tTT and n2 namely -y.^r 



] imply the helpful equation 

^2- - 4l- It is worth 
to note that in classical LVM Eq. [14] similar relation 
exists: = = ^, so we conclude that when n\ ^ 1, 

ri2 ^ 1 the results of statistical description completely 
coincide with pure dynamical consideration. But in the 
case of small numbers ni, ni the difference between them 
may be essential. To demonstrate this distinction and 
also to compare the approach proposed in the present 
paper with usual Markovian description of such systems 
proposed in well known article of Nicolis and Prigogine 



@ and expanded in their later book @ it is appropriate 
to consider the truncated case of LVM, with Ai = A2 = 0. 
It is obvious that total number of individuals N= n\ +?i2 
in this model will be constant .This fact greatly simplifies 
finding and analysis of solutions Eq. [TH]which in this case 
takes the form: 



dG_ 
'dt 



{v 



dudv 



(vG) 



(21) 



It is easy to see that Eq. [21] has solutions for any inte- 
ger N in the form of homogeneous polynomial in u and 
V of degree N , namely Gat {u,v,t) ~ J2^k {t)u''v^~'' 

k 

where coefficients Ak satisfy the normalization condition 
J2^k = 1- Thereby Eq. [5T] is reduced to the linear 
fc 
system of equations of N-|-lorder for coefficients Ak of 

the form —r^ = LkmAm, where matrix elements Lkm 
can directly be found from Eq. [5T]for any N. For exam- 

-2 the matrix Lhrn has the 



pie in simplest case when N= 

/-2 0^ 
form 2-2 

\ 2 Oy 

Let us consider now this case more detail. Let 
G2 (u,v,t) — av? + buy + cv^ is generating function of 
the model. Then Eq. [21] implies the next system of 
equations for evolution of coefficients a, &, c. 



da 
'dt 



-2a, 



db 
di 



2a -2b, 



dc 
di 



2b. 



(22) 



Together with normalization condition a + b + c = 1 sys- 
tem Eq. [22] allows one to give statistical description of 
the model at any time. In particular Eq. 1221 implies that 
when t tends to infinity riT tends to zero and nT tends 
to 2. This completely agrees with solutions of dynamical 
equations in this case. From the other hand let us con- 
sider now the equation for generating function of this 
model obtained by Nicolis and Prigogine (see eq .10.66 
in their book [6|) which in our notation has the form 



dG 
'dt 



u)v 



d^G 
dudv 



(23) 



In the case N — 2 Eq. [251 implies the system of equations 



for coefficients of expansion G2 = au^ 
differs from [22] namely 



buv 



da 
'di 



db 
'di 



dc 
'di 



which 



(24) 



Eq. ([24| imply that when t tends to infinity both quan- 
tities rii, fL2 tend to nonzero values depending on initial 
conditions what obviously does not agree with dynami- 
cal equations. But when A^ ^ 1 it is easy to see that 
asymptotic behaviour solutions obtained from equations 
Eq. dSl]) and Eq. ^^ are virtually identical. Thus 
in this example we see the benefits of approach proposed 
over the standard methods which do not allow to consider 
explicitly the discreteness of variables of the problem. 

Finally in the last part of the paper using one concrete 
model of the population dynamics we want to demon- 
strate that approach proposed let one to give statistical 
description of the systems for which dynamical descrip- 
tion in the framework of mean field theory looks as over- 
simplified. For this purpose we consider two "compet- 
ing" kins of cannibals eating each other so that voracity 
of individuals in both kins is assumed to be distinctive. 
As one knows cannibalism is widespread in living na- 
ture and plays important role in evolution processes [7|. 
Besides many species possess special mechanisms which 
let them to recognize relatives and to avoid of their eat- 
ing Q. Let us assume that mutual eating is the major 
factor of changes of the number of individuals in both 
kins. Then evolution of the number of individuals niand 
?i2 in such model can be represented of simple system of 
equation of the form 



drii 
IT 
dn2 
~df 



anin2 — bnin2, 
—anin2 -I- bnin2. 



(25) 



It is easy to see that total number of individuals in this 
model N= ni -I- n2 conserves. But dynamical descrip- 
tion of the system with the help of Eq. ([25|) seems to 
be oversimplified. In particular it implies that for any N 
when a > b and t tends to infinity rii tends to N , and 
n2 tends to zero. Now following the spirit of our method 
we will describe this system by the help of two operators 

i?i = AioiaJ and R2 = X2afa2- Show that such statis- 
tical version is completely consistent. Actually acting as 
in previous examples we can write the master equation 
for the density matrix of the system as 



dp 
'dt 



Rip^Ri + i?2P, ^2 



(26) 



If again we are interesting in by solutions of the Eq. 
(^5)1 of the form p = X) l^i'^a) P„„„2 ("i'^2| then for 

nin2 

the generating function of the problem G{u,v,t) — 
S Pnin2'""^'^"^ ^^ obtain the equation 



fir H'^ r/2 

a{u~v)-—{uG) + b{v-u)-—{vG), (27) 
at cluav duav 



where a = 2\\ and b = 2X\. The equation Eq. ((27)) 
implies that generating function of stationary state of 
this model in the case when total number of individuals 
is equal to N can be represented in the form 



Gst{u,v)^C^ 



N 



(bv) 



N+1 



(au) 



N+l 



bv 



C 



N 



{bvf + {bvf~^ (au) + , 



(28) 



r 



the 

, -1 



normalization factor, G, 



N 



where Gn is 

(6^ + ab'^~^ + ...) ^ Having in hands expression Eq. 
(PS)) for the generation function one can find all statisti- 
cal characteristics of the model for any N. In particular 
for N=2 for the average values of individuals in both kins 
we obtain nT = J°:^t°'^^ , and rT7 = jft?^ . Now we 



b'^+ab+a^ ' 



h2+af,+a2 • 



want to show that when N tends to infinity the proper- 
ties of our statistical model will be similar to ones of the 
dynamical model Eq. ([25]) . Let us assume that b > a and 
let X — ^ , {k < 1) . Let us calculate now the ratio 

= in this model. Using Eq. (^5]) we obtain 



W _ Na^ + (iV - 1) a^-^b + 
^^ Nb^ + {N-l)b^-^a + 






(29) 



Expression Eq. 



can be represented in the next convenient form 



7T4 



iV- 



,2 dLnfNJx:) 
"^ a 

OK 



^ (9Lk/jv(>c) 



(30) 



where /at (h) 

(JV-l)>f"-JV> " 



-1 



-. It is easy to see that ^ g^" — 

(i->c)(i->c'^) • Obviously this expression tends to 

j^ when K < 1 and N tends to infinity .Thus Eq. ([50]) 

imphes that ratio = tends to zero but this case is reahzed 
only for infinite population. From the other hand for any 
finite population of cannibals all its statistical properties 
can be well described by the approach proposed in our pa- 
per. In conclusion let us briefly summing up our consid- 
eration. We propose the statistical method of describing 
of various systems in physics, chemistry,ecology whose 
states can be represented by integers. Although the 
method proposed is based on quantum Lindblad equa- 
tion nevertheless we have shown that it can be success- 
fully used also for description of classical open systems 



with integer variables (at least in the cases when "Hamil- 
tonian" of the open system is equal to zero) .All examples 
considered in present paper demonstrate that approach 
proposed results in consistent conclusions and allows one 
in several cases to eliminate essential inaccuracies of pre- 
ceding considerations. 
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